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DEVELOPMENT OF A FRACTIONAL-STEP METHOD FOR THE UNSTEADY 
INCOMPRESSIBLE NAVIER-STOKES EQUATIONS IN GENERALIZED 

COORDINATE SYSTEMS 


Moshe Rosenfeld,* Dochan Kwak, and Marcel Vinokurt 
Ames Research Center 

SUMMARY 

A fractional step method is developed for solving the time-dependent three-dimensional incompress- 
ible Navier-Stokes equations in generalized coordinate systems. The primitive variable formulation uses 
the pressure, defined at the center of the computational cell, and the volume fluxes across the faces 
of the cells as the dependent variables, instead of the Cartesian components of the velocity. This 
choice is equivalent to using the contravariant velocity components in a staggered grid multiplied by 
the volume of the computational cell. The governing equations are discretized by finite volumes using 
a staggered mesh system. The solution of the continuity equation is decoupled from the momentum 
equations by a fractional step method which enforces mass conservation by solving a Poisson equation. 
This procedure, combined with the consistent approximations of the geometric quantities, is done to 
satisfy the discretized mass conservation equation to machine accuracy, as well as to gain the favorable 
convergence properties of the Poisson solver. The momentum equations are solved by an approximate 
factorization method, and a novel ZEBRA scheme with four-color ordering is devised for the efficient 
solution of the Poisson equation. Several two- and three-dimensional laminar test cases are computed 
and compared with other numerical and experimental results to validate the solution method. Good 
agreement is obtained in all cases. 


1. INTRODUCTION 

In simulating viscous incompressible flows, much effort has been directed to obtaining steady state 
solutions (refs. 1-3), partly because of limited computer resources. However, numerous interesting and 
important fluid flow phenomena are essentially time dependent, i.e., flow separation, vortex shedding, 
turbulence, turbomachinery flow and flow in biological systems. Numerous time-dependent solutions 
of viscous incompressible flows have been reported for two-dimensional cases (refs. 4-6). However, 
to obtain time-dependent solutions of realistic three-dimensional problems within a reasonable time, it 
is essential to have a computationally efficient flow solver. The need for such a general solver has 
been recognized, but only recently, with the advent of new-generation super-computers, has it become 
an achievable goal. To date, only a limited number of unsteady, three-dimensional solutions of the 
incompressible Navier-Stokes equations in generalized coordinate systems have been reported (ref. 7). 
Therefore, the purpose of the present study is to develop and validate an accurate, unsteady, viscous, 
incompressible flow solver for arbitrary geometries. 
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Direct solution of the discrete equations resulting from the incompressible Navier-Stokes equations 
is still impractical for general three-dimensional cases. An iterative solution is difficult because of 
the lack of a pressure time-derivative term in the continuity equation. Several approaches have been 
suggested to overcome this problem. In the pseudo-compressibility method (refs. 1 and 6-8) a fictitious 
time derivative of the pressure is added to the continuity equation sO that the modified set of equations 
can be solved implicitly by marching in time. When a steady solution is reached, the original equations 
are recovered. Time accuracy can be achieved by using a dual time-stepping (ref. 8) or a pseudo time 
iterative procedure (refs. 6 and 7), in which the solution at each time step is obtained by sub-iterations 
on a set of equations similar to what is devised for the steady pseudo-compressibility method. 

Another approach used in many incompressible flow computations is the fractional step (projection) 
method with its numerous variants (refs. 9-14). The solution is advanced one time step in two (or more) 
stages. Usually, in the first stage the momentum equations are solved for an approximate velocity field 
which is not generally divergence-free. In the second stage, the pressure and the velocity fields are 
corrected to satisfy the continuity equation. This step leads to a Poisson equation with Neumann-type 
boundary conditions. An efficient solution of the Poisson equation is critical since it may consume a 
substantial portion of the total computing time. The present procedure is developed based on a fractional 
step approach. The variable definitions and the staggered mesh arrangement are chosen to facilitate the 
development of an efficient Poisson solver for curvilinear coordinates. A novel and efficient ZEBRA 
scheme with four-color ordering is devised for the nonorthogonal Poisson solver. Using this scheme, a 
multigrid acceleration procedure can be readily incorporated. 

An important aspect of achieving accuracy in arbitrary domains is related to the method of dis- 
cretization. For example, certain geometric identities have to be satisfied accurately in the discrete sense 
as well as in the continuous domain. In this respect, the finite-volume approach can more easily yield 
accurate and conservative approximations than methods based on finite -differences (ref. 15). Therefore, 
a finite-volume discretization method in a staggered grid has been used with a consistent approximation 
of geometric quantities in generalized curvilinear coordinate systems. 

Another important aspect related to the efficiency and the accuracy of the solution method is the 
choice of the dependent variables. The difficulties associated with the pressure may be eliminated 
by replacing the pressure with the vorticity vector. However, the order of the governing equations 
is increased with adverse effects on the computational time requirement, and the specification of the 
vorticity boundary conditions is not straightforward. The primitive variables are the natural dependent 
variables for general three-dimensional flow calculations. Most existing primitive-variable solution 
methods use the Cartesian velocity components to describe the velocity vector. Shyy et al. (ref. 3) used 
contravariant-type velocity components, but only in the second step of a fractional step method (the 
correction stage). In incompressible flows, most of the solution methods use a staggered mesh in order 
to prevent the problems associated with the “checkerboard” behavior of the pressure field (ref. 13). The 
extension of this approach to generalized coordinate systems is not straightforward, since the Cartesian 
components of the velocity cannot be uniquely related to particular coordinate lines. 

In the present study, the volume fluxes across each face of the computational cells are chosen as 
dependent variables, instead of the Cartesian components of the velocity. These variables correspond to 
the contravariant velocity components in a staggered grid multiplied by the volume of the computational 
cell, and they make possible a simple extension of the staggered grid approach to generalized curvilinear 
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coordinate systems. Using this choice, the discretized mass conservation equation can easily be satisfied, 
with favorable effects on the convergence properties of the Poisson equation solver. 

The governing equations in an integral form are given in Section 2. The finite-volume discretization 
procedure and the numerical solution method are elaborated in Sections 3 and 4, respectively. Com- 
puted results of several validation cases are compared with other numerical and experimental results in 
Section 5, followed by a summary in Section 6. 

2. FORMULATION 

The equations governing the flow of constant density isothermal incompressible fluids in a fixed 
control volume V with face S, are the conservation of mass 

<j> s dS - u = 0 ( 1 ) 

and the conservation of momentum 

ltfv udV = i dsT m 

where u is the velocity vector, dV is a volume element and dS is the area element vector. The tensor 
T is given by 

f — -uu - P I + v [Vu + ( Vu) T ] ( 3 ) 

for Newtonian fluids, where / is the identity tensor, Vu is the gradient of u, and (')^ is the transpose 
operator. The pressure is P and u is the effective kinematic viscosity. The viscous terms are written in 
a form easily replaceable by a variable kinematic viscosity coefficient (turbulence models) or by certain 
non-Newtonian flow models. 

By using the divergence theorem, equations (1) and (2) can be written in an equivalent differential 
form. The integral formulation is preferred in the present work since a finite-volume discretization 
method is employed. 


3. DISCRETIZATION 
Geometric Quantities 

A general nonorthogonal coordinate system (£, rj, Q is defined (discretely) by 

r = r(£,r?,C) (4) 

where r = ( x,y,z) T is the Cartesian coordinate system. In the present work, only fixed coordinate 
systems are considered, although the extension to moving coordinate systems is possible (ref. 16). The 
computational domain (£, r ], £) is divided into uniform primary cells with mesh size A£ = A r\ = A£ = 
1 , and the center of each primary cell corresponds to the indices i, j, k. In the finite-volume discretization 
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procedure, the integral governing equations are approximated over the computational volumes in the 
physical space. 


The face l of a cell is given by the vector quantity (fig. 1) 


d(l+ 1) d(l + 2) 

The computational coordinates l = £,r) or C are in cyclic order and x is the cross product operator. The 
vector quantity has the magnitude of the area of the face and a direction normal to it. The differential 

analog is S ; = —V/ , where J is the Jacobian of the inverse of the transformation (4) and V/ is the 
contravariant base vector. 

As Vinokur (ref. 15) pointed out, an accurate discretization should satisfy certain geometric con- 
servation laws. The condition that a cell is closed (a special case of (1», 

dS = 0 (6a) 

should be satisfied exactly in the discrete form as well: 

E s ' = ° (6b) 

faces 

where the summation (with the proper signs) is over all the faces of the computational cell. Equation 
(6b) may be satisfied by computing S l at the center of each face by a proper approximation of — in 

t/i 

eq. (5). For example, the area vector S^, which is related to the point ( i + may be computed 

by using the second-order approximation 

(7) 

~ ?{ r j- T j~2,k~2 + r j+2,k+\ r j+^,k-j)i+^ 

It is assumed that the coordinates of the cell vertices are given. 

In order to ensure that the volumes of all the computational cells sum to the total volume, the 
volume of each cell is computed by dividing the cell into three pyramids having in common the main 
diagonal and one vertex of the cell (ref. 15): 


, H~ S T> , + s? 


r i+3ij+3>fc+2 — 2>j~ 2ifc~ 2 


The volume V of the cell is the inverse of the Jacobian: V t j) c = 1 / Jj 
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Mass Conservation 


Discretization of the mass conservation equation (1) for a computational cell which coincides with 
the primary cell yields 


(S { ' “),+ ■ - (S« • u) H + (S’ . u) , - (S’ ■ «) ..,+ 

(9) 

(S C • »)*+! - (S C • u) fc _ , = Y faces S l ■ U = o 

where Y implies summation (with the proper signs) over all the faces of a computational cell. Note that 
throughout the present paper the indices (i,j, k ) are omitted for simplicity whenever possible. Each term 
on the left-hand side of (9) approximates the volume flux over the corresponding face. Equation (9) 
states that the net mass flux (divided by the constant density) over each cell is zero, as no mass is 
generated within the cell. A discretized mass conservation equation, which is identical in form to the 
Cartesian case, can be derived if the volume fluxes over the faces of the computational cells are chosen 
as the unknowns instead of the Cartesian velocity components. Let 


U* = • u 

U 11 =S' n -u (10) 

(/ c = S c • u 


where L/ 77 and are the volume fluxes over the £, rj and £ faces of a primary cell. Then, the 
continuity equation in any coordinate system takes the simple form 



-vf_ i + vl 1 -ul , + c£ + , - ui , 






k-i 


= D iv (U l ) = 0 


( 11 ) 


where U l = ({/£, U 71 , U *>). The summation operator D{ v is a discrete divergence-like operator (the di- 
vergence operator itself is —D( v ). In tensor algebra nomenclature, U l are the contravariant components 

of the velocity vector (in a staggered arrangement) multiplied by the volume V of the corresponding 
computational cell. 


Accumulated experience with fractional step solution methods of unsteady flow shows the impor- 
tance of exactly satisfying the discrete mass conservation equation (ref. 10). Therefore, the simple form 
of (11), which can be satisfied to round-off errors in any generalized coordinate system, suggests that the 
volume fluxes are the natural dependent variables in the context of fractional step methods. This choice 
complicates somewhat the discretization of the momentum equations, but is important for obtaining a 
divergence-free velocity field in generalized coordinate systems. 

Summing (11) over all cells yields the global mass conservation equation 

E^=0 (12) 

B 

where Yb implies summation (with the proper signs) over all the boundaries of the computational 
domain. If velocity components are specified as the boundary conditions over all the boundaries, 
consistency requires the exact satisfaction of the discrete equation (12). 
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Momentum Conservation 


Spatial discretization of the integral momentum conservation equation (2) applied on a cell with 
constant volume V yields 

du , i 

V— =ys l -T = F (13) 

a ^ 

The momentum equations should be reformulated so that the volume fluxes will be the dependent 
variables. This can be done by a scalar multiplication of the vector momentum equations with the 
corresponding face area vector S l . Hence, the momentum equations are projected onto the direction of 
the corresponding face area vector S* (and multiplied by its magnitude). This projection may have a 
favorable numerical effect provided the coordinate system is approximately aligned with the flow field. 
The discretization of the momentum equations for each U * is performed on a different computational 
cell. Each cell has the formal size of A£ x A 77 x AC in the computational space, but the centers are 
located at (i + j,j,k), ( i,j + 2 ,A;), anc * (*>.?»& + 2 ) for the U^U 71 , and momentum equations, 
respectively. This choice is in accordance with the conventional staggered grid discretization practice. 


The derivation of the ^-momentum equation will be described in this section. The other two 
momentum equations can be obtained by cyclic permutation. 

The control volume used for the discretization of the [/^-momentum equation is shown as the 
shaded region in figure 2. The center of the cell is marked by the indices ( i + j,k). Dotting 

with (13) yields 


i+i 


s* 


t-H 


V 


du\ 

Tt) l 


= v. 


i+2 


2 + 


dl/t , 

dt 


= S* , 

x+2 




( 14 ) 


for a fixed grid. The term F ;4 j is the total flux through the computational cell for the ^-momentum 
equation. The evaluation of 

scalar version since the flux vector should be computed on each face. 


equation. The evaluation of F 1 requires three times as much computational work as the standard 

(> 1 J 


The static pressure P and the effective kinematic viscosity u are defined at the center ( i,j,k ) 
of each primary cell. The definition of the discrete pressure at this point is important for obtaining 
the standard second-order approximation of the Poisson equation, which is derived in the process of 
the fractional step solution of the discrete equations (see Section 4). The present definition of the 
volume fluxes and the mid-cell location of the pressure are equivalent to the staggered formulation of 
conventional fractional step solution methods. 

In the following subsections, the temporal and the spatial discretization of the (/^-momentum 
equation will be elaborated separately. 

Temporal discretization- The (7^ -momentum equation (14) is rewritten in operator form as 


dU C , 


(15) 
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where L £ includes the right-hand-side terms of equation (14). For numerical reasons, the operator 
L(: must be split into several parts, which should be treated differently in the process of the time- 
discretization and the numerical solution 

h = C (P) + C(,e(U l ) + R((P) + D ( (U‘) + D te (U l ) (16) 

The convection terms are split into the nonlinear operators C^ e {U l ) and C^(U l ) which include the 
explicit and implicit terms, respectively. The linear operators D^(U l ) and D^ e (U l ) are the implicit 
and explicit parts of the diffusion terms, respectively, while the operator R^(P) includes the pressure 
terms. The approximate factorization solution method (see section 4) and the numerical solution of the 
momentum equations require explicit treatment of certain convection and diffusion terms, which do not 
fit a scalar tri-diagonal pattern of the coefficient matrix of the approximate factorization steps. 

A general second-order-accurate scheme, with explicit approximation of certain convection and 
diffusion terms, is given by the three-time-level scheme 

(1 + e)(^) n+1 - (1 + 2e)((7^) n + e^)"- 1 _ 

i+ h At ~ 

(| + e)ty - (j + e)L n f' +0 r (rt"+‘ - 2if + ft"- 1 ) + (17) 

«c(C"+' -2C? \- C"-' j + 0 d (D"+' -2 Dj? + £>"-’) 

where n is the time level and e. 6 r . 6 C and 9 d are parameters to be selected for a particular scheme. 
The partly explicit approximation of certain convection and diffusion terms may restrict the maximal 
allowable time step for stable solutions in highly non-orthogonal meshes. Yet, this restriction is usually 
acceptable for time-dependent solutions where physical reasons limit the time step anyway. 

One of the simplest choices for the free parameters is 9 d = 9 r — 0.5, 9 C = 0 and e = 0 (Adams- 
Bashforth/Crank-Nicolson scheme). However, an analysis by Beam and Warming (ref. 17) demonstrates 
that the split form of a purely parabolic equation with explicit approximation of the mixed derivatives is 
not stable for this choice of parameters. Stable solutions may be achieved by a proper choice of e, 9 r , 9 d 
and 9 C . Generally, the analysis of a simplified parabolic model problem shows that a necessary condition 

2( 1 + f)^ 

for the stability of the factored scheme is e > — 1/2, and 9 > ^ - - (for 9 d = 9 r = 9 C — 9) (ref. 17). 

A more general scheme can be obtained by choosing A — (AU^) n — a u (AU^) n ~^ and 
A P — A P n — a p AP n ~ l as the unknowns. Note that Ad> n = <f> n+1 — $ n and $ = (U^,P). 
Equation (17) for A U ^ and A P takes the form 

K + , A US = G + fai A U l ) + 9 r R^A P)) (18a) 

where 

E, c =9 c C^(AU l ) + 9 d D^(AU l ) 

G ={jT~e~ <*u) V. + ± (AUt) n ~' + ^((i + e)q - (i + e)L n ^+) (18b) 

x(a p - l)9 r R^(AP) n - 1 + K - l)E^AU^) n ~ l 
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includes all the terms from the previous two time levels and the parameters a u ,ct p are to be selected. 
The choice a u = a p = 0 corresponds to the well known “Delta form” where the unknowns are of 
order At, whereas a u = a p = 1 results in the most accurate scheme for an approximate factorization 
solution method since the unknowns are of order {At) 2 . 

Enhanced stability properties can be obtained by an implicit approximation of the convection terms 
(9 C > 0). The implicit convection terms are linearized in the present study to first-order accuracy in 
time. For example, a general term U ■ U is approximated by ( U ■ U) n+l « U n ■ U n+ . To retain the 
scalar tri-diagonal (or penta-diagonal) form of the approximate factorization steps (see section 4), only 
the convection terms that fit into this pattern are approximated implicitly. 


Spatial discretization- The finite-volume discretization of the ^-momentum equation is performed 
for a cell centered at (i + j, j, k), as is illustrated in figure 2. The time marching scheme (18) requires 
the computation of the operator L ^ and of the operators and which also appear in the implicit 
part of (18a). For the time level n, L ^ is defined by 


L? = s e , -J2S 1 -T n (19) 

1+2 l 

where the summation (with the proper signs) is over all the faces of the computational cell. To compute 
the fluxes over each face, the velocity vector should be computed from the volume-flux unknowns by 
the identity 

u = 4- S V U V + S C U C = S m U m (20) 

From (10) and (20) one obtains 

U l = S l ■ u = S ; ■ SmU m (21a) 

The invariance of the velocity vector requires 

S*-S m = 6 l m (21b) 

where <5^ is the Kronecker delta and S m is the inverse base of S^. It has the differential analogue 

(9r ^ 

S m = J — — so that in terms of tensor algebra, S m is the covariant base vector multiplied by the 

dm 

Jacobian J, while S^ is the contravariant base vector multiplied by 1 j J . A uniform velocity field can 
be numerically preserved if the base S m is computed at each point from the relation 


gm+l x gm+2 

Sm x gm+2) (22) 

which satisfies (21b) identically. Here m is the cyclic permutation of (£, rj, Q. 

The £-face center of the L^-momentum cell is at the point (i, j, k) (fig. 2). The flux over this face 
is computed from 


(S* • T) iJ>k = (-UtySt -Stp + St-is (Vu + (Vu) T )) (23) 

V ' i,j,k 


8 


I 



The conservative form of the velocity vector gradient is 


j>^d Su 


Applying (24) to the computation of yields 


VUj i k = l U • , 1 — ,U. l 

l 'J' K V\ l+ 1 l- 2 


+S 77 , u. i - S v ,u ! +$l ,iU. ,1 -Sj* ,u , 

J + 2 3+1 J - 2 -? - 2 k +l fc+ 2 fc -2 2 
= 77^5, iS, ■ , i u\ i — sf jS, . _\U. i 




+s U s ^ u k 


u 1 ., - S 77 ,S. 

J + 2 ^3 7 3 7 


+sf 1 S.^,^Li-Sj ,S.._ iC/f 1 


The fluxes at the point (i, j, k ) are computed by a second-order-accurate averaging and therefore 
the scheme is equivalent to second-order central differences. The details of the computation of the 
geometric quantities, the volume fluxes and the pressure and kinematic viscosity at the various points 
of the staggered grid, are given in Appendices A and B. 


The flux of the r^-face is approximated at {i + \,j — \,k) similarly by 


< s ” • fy+ij-i* 


= (-U 1] U l Si - S V P + S 11 -V (Vu + (Vu) T )] 




where the gradient is computed from 


^ U »+zJ“2>* V. \ . 1 


s*s iU l ), 


1 - {S^SiU 1 ). 1 
i+ Ij — 2 3~1 


+(S , S ,(/'), + , - (S*s,c/‘) i+iJ _ 


+(S<S,U‘) i 


i+\,j—?,k+l 


The flux of the C-face is approximated at (i + j, j, k — j) by 


(S<S ,U l ) 




(S C • f )i + Uk-i = (~^U% -S<P + S<-v (Vu + (Vu) T )) j 


1 „• u i 
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where the gradient is computed from 


h + L k -i (&SiU l )i 




(S«S 


+(S”S ^'), + . - (SIS ,u‘) i+hH , k _, 


+(S<S,(7') i+ j 


(S<S ,u l 


i 4 - 7 , fc — 1 


Computation of - summary: The computation of the operator L ^ (see definition in (19)) is 
summarized below. 

For each face l : 


(a) compute the tensor Vu. 

(b) compute the transpose (Vu) T . 

(c) compute the viscous flux • v (Vu + (Vu) T ). 

(d) add the inviscid flux - U l U m S m - S l P 

(e) store result in G*. 

(f) dot with S e , : N l = S ? , • G l 

i+\ i+3 


For each cell compute 



+ N V j . , 

J+? 


N'\ , . , + N C . 



1 

1 


Definition of the operator R\\ The gradient-like operator R[ is given by R\ = — V m+ 1 S'- , • VP 

(no summation is implied), where VP is the gradient of P. It is proportional to the projection of the 
pressure gradient on the direction normal to the l cell face. Appendix C gives the terms of Rp 

Definition of the operators D[ and C\\ The operators D/ and Q include the implicit terms that fit 
the pattern required by the approximate factorization method used for solving the momentum equations 
(see section 4). For nearly orthogonal coordinate systems, the most significant diffusion and convection 
terms are included in the implicit part D\ and C\. Yet, this choice requires the solution of a block- 
tridiagonal system because the momentum equations are coupled. The equations can be decoupled by 
calculating implicitly the terms that involve only the unknown of that momentum equation (e.g., the 
terms of U ^ in the momentum equation). The explicit computation of the other terms may reduce 
stability (ref. 17), but the simplicity obtained may be valuable for time-dependent solutions. It should 
be recalled that precautions were already taken to enhance stability of such schemes (see the section of 
Temporal Discretization). The definition of the D £ and operators is given in Appendices D and E, 
respectively. 
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In certain high Reynolds number cases, artificial diffusion should be added to smooth out the 
solution. In the present work, fourth-order numerical diffusion is used by adding it directly to the 
fluxes at the faces of the computational faces. Lower-order numerical diffusion is not used to maintain 
the second-order accuracy of the scheme and to minimize the effect of the numerical diffusion. The 
fourth-order numerical diffusion may be added explicitly or implicitly. 


Geometrical considerations- The approximation of the fluxes at the computational cell faces re- 
quires the determination of the geometric parameters and S i at numerous points for each computational 
cell. If uniform velocity is to be preserved, S; should be evaluated at each point from (22) at the cost of 
more than 2000 additional operations for each cell. On the other hand, the available memory of present 
computers still imposes certain limitations on the number of geometric parameters that can be stored. In 
the present work, a compromise between CPU time and storage requirements has been made by storing 


only the quantities S . 


The base is computed from (5) while S i is computed from (22). 


*~ 2 ’ 


3 - 


1 > 
2 








k _ i and Vi j f. for all computational cells. 


The geometric quantities at the other points are computed by simple averaging, see Appendix A. 
It can be shown (ref. 15) that this averaging satisfies the geometric conservation laws (closed cells and 
invariance of the total volume) even though a uniform free-stream velocity is not fully preserved since 
(22) is not satisfied at all points. In most cases the approximation of the S; terms by simple averaging 
will not lead to a serious degradation of the results. This approximation affects only the computation 
of the momentum equation fluxes (the continuity equation is solved to machine accuracy). Near solid 
boundaries where large gradients usually occur, this approximation introduces only small errors, since 
a fine (and nearly orthogonal) mesh can be used. Far from a solid boundary, the velocity gradients are 
small (especially for external flows) so the errors introduced by this approximation are small. Moreover, 
in these regions the flow is usually potential and therefore the continuity equation, which is not affected 
by this approximation at all, plays the major role. The accuracy of the approximation of the momentum 
equation is a less important factor in these regions. 


4. SOLUTION METHOD 


The Fractional Step Method 

The simultaneous solution of the large number of resulting discrete equations is very costly, espe- 
cially for three-dimensional cases. An efficient approximate solution can be obtained by decoupling the 
solution of the momentum equations from the solution of the continuity equation by a fractional step 
method. The basic fractional step (or projection) method was proposed by Chorin (ref. 9). The MAC 
method proposed earlier by Harlow and Welsh (ref. 14) is actually a variant of that method. The two 
methods are identical as long as the boundary conditions are not considered (see ref. 18) for a detailed 
discussion of the relationship between the numerous variants of the fractional step method). 

In the present implementation, first the momentum equations are solved for an approximate A U l 
by dropping A P ) from (18), so that the pressure gradient is taken from the previous time-level 

( v -n + p--rh E ‘) A& ‘ = G < 30a > 
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where I is the identity operator. The resulting flow field does not generally satisfy the mass conservation 
equation. In the second step A U l is updated by 

v m+i _ (a U‘ - A U l ) = (30b) 

so that the continuity equation (cast in terms of A U l ) will be satisfied at the level n + 1 

D iv [( U l ) n + a u {AU l ) n ~ l + A U 1 } = 0 (30c) 

where m = i,j or k for Z = £,rj, or (, respectively, and Dj v is the summation operator defined by (11). 
The variable (p is a scalar to be defined later. 

Equations (30b) and (30c) are combined into a single discrete Poisson equation, which is easier to 
solve than the coupled set of the equations (30b-c), by applying the operator D{ v on (30b) 


1 + £ 
0 r At 


D 


IV 


(U l ) n + a u {AU l ) 


l\n — 1 


+ A U l 


~~ Div 


m A 

V m + J 


(30d) 


since Di v ((U l ) n+i ) = 0. Note that the discrete Laplacian operator is given by ——D lv {~). Finally, 
the variables at the new time level n + 1 are computed from 


(J7 / ) n+1 = (U l ) n + a u (A U l ) n ~ x + A U l 

A P = 0 (30e) 

pn+l = pn + QpA pn-l +A p 


Substituting A U 1 from (30b) into (30a) does not recover the original discrete momentum equa- 
tion (18) if the substitution A P = <p is used, since the computation of the viscous terms in (30a) is 

based on A U l rather than on A U l . In the case of the fully explicit approximation of the convection 
terms, it is shown in Appendix F that (U l ) n+l is the exact solution of the discrete equations, but 
A P = 4> + 0(uAt) and therefore _pk+* j s not ^ exact discrete solution. For a Cartesian case with 

1 i/Ai 9 

6=0 and 9(i = Kim and Moin (ref. 11) have derived the exact relation A P = 4> —V (f), where 

2 2 


V 2 is the Laplacian operator. However, such a simple relationship cannot be found in generalized co- 
ordinate systems with the present splitting of the diffusion terms (which does not satisfy certain vector 
identities). Braza et al. (ref. 19) have found insignificant differences between solutions which use 
the approximation A P = <f>, and solutions which compute the discrete pressure exactly (in the case 
of vortex shedding over a circular cylinder), although the approximate solution might be less stable. 
Because the difference between 4> and A P is proportional to u, the direct substitution A P = 0 is 
reasonable for high-Reynolds-number flows. This approximation degrades the second-order temporal 
accuracy of the pressure to first-order accuracy since this substitution is of the order 0(At ). However, 
as can be seen from (30b), the velocity solution remains second-order-accurate in time. 
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In the case of the partially explicit scheme ( 9 C > 0), the relation between the pressure and the 
function is more complicated. In this case, the velocity field is also first-order accurate in time because 
of the first-order linearization. 

The first step (30a) is a consistent approximation of the momentum equations, i.e., as (At, A l) — > 
0, (U l ) n+l -* U l , where U l is the solution of the continuous problem. Therefore, the physical 
boundary conditions may be specified from the n + 1 time level. In some fractional step methods the 
sub-steps are not consistent and special boundary conditions should be devised for the intermediate steps 
(ref. 11). 

The drawbacks of the present fractional step method are the increased storage required for the three- 
time-level method and the restriction on the maximum Courant number (CFL) in the case of the explicit 
treatment of the convection terms. Even in the case of the implicit approximation of the convection 
terms, the scheme is not unconditionally stable in the generalized non-orthogonal case because only a 
part of the convection terms is implicitly approximated. The CFL number restriction is tolerable for 
many time -dependent solutions where the physics dictates a time step of the order CFL = 1. 

Solution of the momentum equations- The modified discrete momentum equations (30a) are 
solved by an approximate factorization method. The implicit operator E[ is chosen such that the 
momentum equations are decoupled from one another, and the operator can be decomposed into the 
three parts 

Ei = E l4 + E iiT] + E u (31) 

where Ei ^, Ei rj and E[^ are tri- or penta-diagonal operators in the respective directions. The approx- 
imate solution of each A U l — (A U ^ or A f/ 7 ? or A [/*>) momentum equation is given by 

(^+i 1 - A ' u ‘ = V A (32) 

where no summation is implied on l. The splitting error of this approximate factorization scheme is of 
the order (Af) 3 for the general scheme and (At) 4 for the special case a u = 1. 

The explicit treatment of some of the terms might degrade the stability of the method, especially 
for very skewed meshes. In the present method, special measures have been taken to alleviate these 
problems. The three-time-level discretization method improves the stability properties over the more 
commonly used two-time-level methods, especially when mixed derivatives exist. The three-time- 
level scheme also maintains in the case of linear equations second-order-accuracy in time, even for 
formulations which approximate certain terms explicitly (ref. 17). 
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Solution of the Poisson Equation 


The discrete Poisson equation (30d) can be rewritten as 

L{4>) = / (33a) 

where the source term / is given by 

; (( U‘r + Q„( A u‘Y‘- 1 + A &) (33b) 

If the equations (30b) and (30c) are to be satisfied exactly, the Laplacian-like operator L should be 
computed discretely from the identity 


m = D iv (33c) 

The discrete divergence-like operator Di v was defined in (11). The discrete gradient-like operator 
Rl is given in Appendix C. The resultant form for the discrete Laplacian-like operator L is given in 
Appendix G. 

The normal-derivative (Neumann) type boundary condition for each boundary l is computed from 
(30b), 

( A u ‘ - A 6 ‘) (33d) 
At periodic boundaries, periodicity is used. If Dirichlet-type boundary conditions are specified for the 
velocities, condition (33d) is homogeneous. The compatibility condition necessary for the existence of 
a solution to a Poisson equation with derivative boundary conditions over all the boundaries, f = 0. 
is automatically satisfied because of the imposed global discrete mass conservation (12). Note that 
“numerical” boundary conditions for the pressure are not required since the pressure is defined at the 
center of each cell. 

An efficient solution of the Poisson equation is crucial for the efficiency of the whole solution 
method. For a general nonorthogonal coordinate system, the 19-points star discrete equations pose 
difficult challenges in obtaining fast solutions on vector computers. Many iterative methods cannot be 
efficiently vectorized and suffer from degradation of the convergence rate due to the sharp variation of 
the coefficients. The present solution method uses a novel application of a ZEBRA method with four- 
color ordering to decouple the implicit part of the algebraic equations and enable effective vectorization 
of the Poisson solver, without deteriorating the convergence properties. 

The three-dimensional ZEBRA scheme is an iterative solution scheme which solves implicitly all 
the equations along one coordinate line, say along f , as in the Successive Line Over Relaxation (SLOR) 
method. However, the order in which the lines are processed is not the usual lexicographic order (by 
rows or columns), but a “colored” order, devised so that the implicit solution of a line is decoupled 
from the solution of the other lines that belong to the same color. Most existing applications of the 
ZEBRA scheme use a two-color ordering (“Red-Black” schemes). This ordering is inappropriate for 
nonorthogonal cases. In the present method, the points in the (r), Q plane are classified into four groups 
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and a different color label is given to each group, see figure 3. First all the “black” lines are swept 
in a lexicographic order, then the “red” “blue” and “green” lines, respectively. The implicit solution 
of a line is decoupled from the same color lines; for example, when solving for a “black” line, all 
the neighboring lines are of different color - “red,” “blue,” or “green.” This arrangement enables an 
efficient vectorization of the method, even with the 19-point computational stencil. The implicit solution 
along the £ coordinate line may be exploited to enhance convergence of problems with heavily clustered 
mesh points along that direction. 

A mathematical analysis of the ZEBRA method reveals that although each sweep of a color is 
equivalent to a line-Jacobi iteration, one complete iteration has the better convergence properties of the 
SLOR method. Therefore, the convergence properties of the present Poisson solver are similar to the 
SLOR method, but the required CPU time on a vector computer is significantly reduced because of 
the effective vectorization of the ZEBRA method. Moreover, the ZEBRA scheme has good smoothing 
properties and is nearly an optimal relaxation method for using a multigrid acceleration procedure. 

A detailed study on the properties of the Poisson equation solver and the future implementation of 
a multigrid acceleration procedure will be reported in a separate paper. 

Summary of the solution procedure- The complete solution method of the unsteady incompress- 
ible Navier-Stokes equations is summarized below for a constant time step At: 

(a) Set: 


n = 1 
£ — 0 

At/'_, =AP„_,=0 

ip' = 0 

(b) For l = £,r) and 

(1) Compute (L l ) n 

(2) Compute RHS of (30a), 

+ tr( c i£f-(l + £)£?-' 

+(a p - l)9 r fl { (AP)’>-‘ + ( Q „ - OB^AC/f)— 1 ) ) 


if n = 1 then C\ = 1 else C\ 

V.. , 


3 

2 + 


RHS 


-i^(A u i r~ ] 


(3) solve by approximate-factorization: 



r^-Ei) A U l = RHS 
1 +s V 
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(c) Compute RHS of the Poisson equation: 

/ = - A» ((V‘) n + MAC/')"-' + A £/') 

(d) Solve the Poisson equation: L(<p) = f 

(e) For l — £, ?? and (: 

(1) Compute A U l = A U l + Ri(<t>) 

(2) Compute ( U l ) n+l = (U l ) n + a u (AU l ) n ~ ] + A U l 

(3) Compute P n+1 = P n -(- a p AP n ~ 1 + Q 

(f) For l = £,t) and (: 

Lf-' = Lf 
(A U l ) n ~ x = (A U l ) n 
A P n ~ l = A P n 

n = n + \ 

(g) Repeat (b)-(f) for each time step 


5. RESULTS 

Several representative cases have been solved to validate the solution procedure. The test cases 
include the flow in two- and three-dimensional driven cavities; the unsteady flow over a two-dimensional 
circular cylinder, with and without vortex shedding and the three-dimensional flow in a square duct with 
a 90° bend. Several steady cases are chosen as well since most available experimental and numerical 
results are restricted to this regime. To fully validate unsteady solutions, it is necessary to establish data 
bases for time-dependent cases, in particular for three-dimensional configurations. 

The basic ZEBRA solution method of the Poisson equation is implemented using an over-relaxation 
parameter determined experimentally. No difficulties were experienced in converging the Poisson equa- 
tion to machine accuracy in any of the cases; consequently, a divergence-free velocity field could always 
be obtained. The convergence criterion for the maximal residual was usually set to 10~ 8 , which cor- 
responds to not less than a decrease of five-orders of magnitude in the initial residual. That residual 
equals the divergence of the velocity multiplied by the volume of the computational cell (see (30d)). 

A stability condition of the present solution method cannot be analytically derived. Numerical 
experiments show that the fully explicit convection scheme without artificial diffusion is stable for 
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CFL < f(Re), where CFL is the maximal Courant number defined by the maximum value of 

CFL = (s + + 9 A ‘ = (|c/{| + m + ' u<l) v <34) 

over all the computational cells. The Cartesian velocity components are given by u, v, w, and Ax, Ay, Az 
are the physical mesh sizes in the corresponding directions and V is the volume of the computational 
cell. In the absence of fourth-order smoothing terms, the function f(Re) is a monotonically decreasing 
function of the Reynolds number ( not the cell Reynolds number). Usually, f(Re ) = 0.5—1 results in 
stable solutions for Re < 500. The maximum cell Reynolds number is of the order 10 1 — 10 2 in all the 
computations. By adding fourth-order artificial diffusion terms, the time step can be increased beyond 
this limit without introducing numerical instabilities. 

The partial implicit approximation of the convection and diffusion terms results in a scheme which 
is less sensitive to the CFL number restriction. It was found that CFL < 4 yields a stable solution in 
most cases. If orthogonal grids are used, no CFL number restrition are present. 

The set of the free parameters e, 6 r , 0^, 6 C , a u , and a p can be adjusted to enhance the stability 
of the time-advancing scheme. In the present study no attempt was made to find the optimal values. In 
most of the test runs the following values were used: e = 0, 0 r = 0^ — 9 C = l and a u = a p = 0. 

Some convergence tests have been conducted to confirm the accuracy of the scheme by systematic 
refinement of the mesh size and the time-step. Both non-uniform Cartesian and generalized coordinate 
systems have been tested. The discretized equations are second-order-accurate in both space and time. 
However, because of the substitution A P = <f> (see section 4), the temporal accuracy of the pressure is 
decreased to first-order. The explicit-convection scheme was indeed found to be second-order-accurate in 
space, and the second-order temporal accuracy of the velocity has been verified as well. The convergence 

3 

tests show that the pressure is almost --order-accurate in time, which is better than the theoretical 

prediction. The implicit- convection scheme yields a temporal first-order accurate solution of the velocity 
field as well because of the linearization used in the present study. 

The computer code is written for general three-dimensional problems. Two-dimensional cases are 
solved by using two intervals in the third direction and specifying periodic boundary conditions along 
that direction. The solution method consumes from 0.2 to 0.5 ■ 10 -3 CRAY-YMP CPU seconds/mesh- 
point/time-step, depending on the number of iterations required for the convergence of the Poisson 
solver. A substantial reduction in CPU time is anticipated by implementing a multigrid acceleration of 
the Poisson solver. 

Selected results are presented in the following sections. 

Lid-Driven Cavity Flow 

Two-dimensional driven cavity flow- The two-dimensional driven cavity flow is a standard test 
case for new solution methods of the incompressible Navier-Stokes equations. Although the geometry, 
as well as the boundary conditions are simple, the resulting recirculating flowfield, which has no primary 
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direction may be very complicated. A partial review of previous works can be found in Ghia et al. 
(ref. 21), who also computed that flow over a wide range of Reynolds numbers using very fine grids. 


In the present work, the flowfield in a square driven cavity was solved for several Reynolds numbers 
up to Re — 10 3 , but only the case with Re = 10 3 will be elaborated. Zero velocity initial condition 
was given except at the moving boundary x = 1 where v = 1 was specified ( v is the cartesian velocity 
component in the y direction. The problem was solved separately over two cartesian grids with clustering 
toward the walls. The coarse mesh consisted of 21 x 21 x 3 points and the finer mesh, shown in figure 4, 
consisted of 41 x 41 x 3 points with a minimal spacing of 0.0080 near the walls and a maximal spacing 
of 0.0414 at the center of the cavity. Even this finer grid is coarse relative to the 129 x 129 uniform grid 
(with a mesh-size of 0.008) employed by Ghia et al. (ref. 21) or the 321 x 321 mesh used by Vanka 
(ref. 22). Note the somewhat unusual labeling of the coordinate system (fig. 4) in the present work. 

Figures 5(a) and 5(b) show the stream function contours of the solutions obtained for the meshes 
with 21 x 21 x 3 and 41 x 41 x 3 grid points, respectively, at the non-dimensional time of T = 40. 
Figure 5(c) gives the results obtained by Ghia et al. (ref. 21) (identical contour lines are plotted). The 
solution obtained from the coarse grid is unacceptable, but the solution of the 41 x 41 x 3 grid agrees 
well with the computation of Ghia et al. (ref. 21). Figure 6 compares the vorticity contours obtained for 
the finer mesh solution with the results of Ghia et al. The agreement is good except near the primary 
vortex center. The disagreement is probably caused by the relatively coarse spacing of the present 
non-uniform mesh near the center of the cavity (the mesh size there is about 0.04). 

Table 1 compares several properties of the primary and of the lower left and right secondary vortices. 
The agreement is surprisingly good, bearing in mind the relatively coarse grid that was employed in 
the present calculations and the fact that the present solution has not reached yet the steady-state. The 
second right comer vortex, found in reference 21 was not detected in the present computations because 
of the relatively coarse mesh used. Figure 7(a) compares the present results with the results of Ghia et al. 
(ref. 21) for the distribution of the horizontal component of the velocity along the vertical geometric 
centerline ( y = 0.5), and figure 7(b) compares the distribution of the vertical component of the velocity 
along the horizontal centerline x = 0.5. The convergence of the solution to the results of Ghia et al. 
can be observed with the refinement of the mesh. 

Most computations of the driven cavity flow solve for the steady state only, usually by iterating 
or marching in a fictitious time. The present solution was obtained by a time-accurate method and it 
is interesting to track in time the evolution of the primary and secondary vortices. Figures 8 show the 
instantaneous streamline contours for several non-dimensional times, T, and figure 9 gives the motion 
of the vortices center as a function of time. At the very early stages of the impulsive start of the upper 
sliding wall (at x = 1), a clear primary vortex exists near this wall. At T = 3 a small separation bubble can 
be observed not only at the two lower comers of the cavity, but also near the center of the right vertical 
wall. The flow along that wall is similar to the flow along a flat plate with adverse pressure gradient (the 
velocity is decreasing). The next time shown, T = 4, this separation bubble grows substantially, but is 
distinct reattaches quite at a distance from the lower comer vortex. As time advances, the two vortices 
increase in extent and the mid-wall separation bubble moves downward. Between T = 5 and T = 6 the 
flow field topology changes from a structure of four vortices (a primary vortex, two secondary vortices 
at the lower comers of the cavity and the separation bubble at the right vertical wall) into a structure of 
three vortices. The separation bubble merges with the lower right comer vortex to form a larger single 
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vortex. The separation point along the vertical wall moves downward a small distance and reaches its 
steady state value as soon as T = 14. The reattachment point of the right comer vortex moves to the 
left (downstream) a greater distance but reaches its steady state value at about the same time. 

During the evolution of the right comer vortex, the left comer vortex remains very small and almost 
does not grow. It begins to expand only from T = 14 and reaches its final geometric dimensions at 
T = 22. At that time, the geometric parameters (the centers of the vortices and the separation and the 
reattachment points) of all three vortices attain their approximate steady state values. However, the 
intensity of the vortices, described by the stream-function value at the center, continues to increase. 

The evolution of the large secondary vortex at the right comer pushes left the primary vortex as 
can be seen in figure 9. As long as the primary vortex center is far from the lower left center, an almost 
stationary region persists there and consequently, the left comer vortex is very weak. As the primary 
vortex is intensified and moves toward this comer, substantial positive velocity gradient is built up along 
the lower wall with adverse pressure gradient. The separation point of the left comer vortex moves to 
the right side (upstream) and the geometric dimensions of the vortex increase. 

Two-dimensional driven polar cavity flow- A two-dimensional lid-driven flow in a polar cavity 
was studied by Fuchs and Tillmark (ref. 23), both experimentally and numerically. This geometrically 
more challenging flow is computed in the present study. Figure 10 shows the geometry of the problem 
and the 51 x 51 x 3 mesh used in the present computation. The height AB of the cavity is equal to the 
radius of the boundary AD. A polar coordinate system is used with grid points clustered near the walls. 
Zero velocity is specified on all the boundaries, except on the AD boundary, where a unit tangential 
velocity is given. 

The steady flow was computed for two Reynolds numbers. Re — 60 and Re = 350 (based on the 
height of the wall AB and the tangential velocity of the wall AD), for which experimental and numerical 
results are given in (ref. 23). In the present paper, only the results for the higher Reynolds number will 
be described. 

Figure 1 1 gives the velocity-direction plot (all velocity vectors are plotted with equal magnitude) 
for the Re — 350 case. The resulting flow field is similar to the two-dimensional square cavity flow. 
A main vortex and two secondary vortices are formed in the center and at the two comers B and C, 
respectively. The centers of the vortices agree with a flow visualization obtained by Fuchs and Tillmark 
(ref. 23) within an error of less than 2%. Figure 12 compares the radial and circumferential velocity 
components with the experimental results of Fuchs and Tillmark (ref. 23) along the three radial lines 
9 = — 20°, 0° and 20°. Also shown are their numerical results obtained by a stream-function vorticity 
formulation written in a polar coordinate system. Generally, favorable agreement is obtained among all 
the results. Note in particular the very good agreement in the numerical results, although the present 
numerical results are obtained using about 2.5 times fewer mesh points. The grid points are more 
efficiently distributed in the present study than in the uniform grid employed by Fuchs and Tillmark 
(ref. 23). The small deviation of the numerical results from the experimental results is probably because 
of the three-dimensional effects which were found in the experiment (ref. 23). 

Three-dimensional driven cavity flow- The flowfield in a unit cubic driven cavity was solved in 
several studies, i.e., references 24 and 25. The geometry of the cubic cavity is shown in figure 13(a). 
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The upper wall (x = 1) slides at a constant speed of v = 1 along the y-direction. A uniformly distributed 
grid of 31 x 31 x 31 points was used to solve the flowfield at Re — 100 and Re = 1000. The grid 
is identical with the grid employed by Hwang and Huynh (ref. 24) to facilitate a direct comparison. 
The case with Re = 100 was advanced in time until T = 5 and this solution was used as an initial 
condition for the Re — 1000 case, which was terminated at T = 32. In both cases the steady state was 
not reached, although the changes were confined to small regions in the secondary vortices. 

In figure 13 the distribution of the y-component velocity along the geometric centerline (y = z = 0.5) 
is compared with the numerical solutions of Hwang and Huynh (ref. 24) and Cazalbou et al. (ref. 25). 
Good agreement is obtained with both solutions. 

A more detailed comparison is given in figures 14 and 15 for the case of Re = 100 and in figure 16 
for the case of Re — 1000. In these figures, the direction of the velocity (rather than the velocity vector 
itself, which is very small in the comer regions) is compared with the corresponding results of Hwang 
and Huynh (ref. 24) for several cross sections. Figure 14 compares for Re = 100 the velocity direction 
in the symmetry plane (z = 0.5) and in the plane z = 0.033, next to the side wall. The agreement is 
favorable, except near the lower left comer of the plane z = 0.033 where the velocity is very small in 
magnitude. Figure 15 compares the velocity direction for the sections y = 0.63,0.57,0.50,0.47. Again, 
the agreement is good, except at the vicinity of the two comer regions where the velocity magnitude 
is very small. Although the analysis of the flow will not be attempted in this paper, the reader should 
note the very complicated flow patterns generated even in that simple geometry low Reynolds number 
flow case. 

Figure 16 compares the velocity direction of the higher Reynolds number case, (Re — 1000), for 
the sections z = 0.50, 0.033 and y = 0.57. The agreement is satisfactory in this case too, except in the 
prediction of the secondary vortex extent at the lower right comer for 2 — 0.5. The present solution 
predicts a smaller vortex there. It has been noticed that the geometrical extent of this vortex decreased 
with time and this may raise the possibility that neither the solution of Hwang and Huynh (ref. 24) is 
fully converged in this region. 


Flow Over a Circular Cylinder 

The two-dimensional flow over a circular cylinder is solved as an example of an external flow over 
a bluff body. The resulting flow field strongly depends on the Reynolds number. For 40 > Re > 6 
a steady state exists, with a pair of symmetric separation bubbles on the leeward side. At higher 
Reynolds-numbers the flow field is inherently unsteady and is characterized by cyclic vortex shedding. 

In the present work, the time -dependent flow field over a two-dimensional circular cylinder is 
computed for a range of Reynolds numbers 1000 > Re > 40 (based on the diameter). The symmetrical 
flow is solved for an impulsively started cylinder at Re = 40 and Re = 550 while the asymmetric case 
with vortex shedding is simulated for Re = 100, 200, 550, and 1000. A no-slip condition is given on 
the cylinder and uniform velocity is specified on the far-field boundary. 

Symmetric flow- The symmetric flow for Re = 40 is solved using a nonorthogonal coordinate 
system with 45 x 73 x 3 mesh points in the radial, circumferential, and axial directions, respectively. 
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A non-concentric circular outer boundary with a radius of 20 units (cylinder diameters) is constructed. 
The center of the outer boundary is shifted into the wake region 15 units from the center of the cylinder. 
The case with Re = 550 uses a cylindrical grid with a concentric outer boundary at a distance of 
50 diameters and 81 x 85 x 3 mesh points. In both cases, mesh points are clustered near the cylinder 
and in the wake region. 

The time evolution of the separation length (measured from the rear of the cylinder and normalized 
by the diameter) for Re = 40 is compared with the numerical results of Collins and Dennis (ref. 26) 
and with the experimental results of Coutanceau and Bouard (ref. 27) in figure 17. Good agreement is 
obtained, especially at the initial stages of the flow evolution ( t < 8). The slight discrepancy for t > 8 
is attributed to “wall effects” which exist in the experimental case because of the water tunnel walls 
and in the present numerical solution because of the Dirichlet-type boundary conditions at the outer 
boundary. 

At a higher Reynolds number, stable symmetrical flow exists only for a short period of time. In 
figure 18 the time evolution of the separation length for Re = 550 is compared with the numerical 
solution obtained by Loc (ref. 28) using a two-dimensional stream-function vorticity formulation and 
with the experimental results of Bouard and Coutanceau (ref. 29). The two numerical solutions compare 
favorably, but slightly under estimate the experimental results. The good agreement in the numerical 
results is observed although totally different methods and grids are used. In the present study a sub- 
stantially coarser mesh is employed and a spatially second-order-accurate scheme is used, while Loc 
(ref. 28) uses a fourth-order-accurate scheme. 

Asymmetric flow- At a Reynolds number higher than 40 any perturbation excites an unsteady flow 
and eventually a periodic vortex shedding is established generating the well-known von Karman vortex 
street. This phenomenon has been addressed in several previous numerical and experimental works; see 
references 5, 18, 19, and 27 for a more comprehensive review. 

In the present study, the laminar vortex shedding over a circular cylinder is simulated at Reynolds 
numbers of 100, 200, 550 and 1000. The analysis of the flow field properties and the explanation of the 
dynamical phenomena associated with it are beyond the scope of the present paper. The main concern 
here is in validating the numerical method and studying its properties and capabilities. 

A cylindrical coordinate system is used with 81 x 85 x 3 mesh points in the radial, circumferential 
and axial directions, respectively. The concentric outer boundary is at a distance of 50 diameters from 
the center of the cylinder. Mesh points are clustered near the body and in the wake region. The minimal 
radial spacing near the cylinder is 0.014, while the maximal radial spacing at the outer boundary is 2.15. 
The ratio between the maximal and the minimal Jacobian is of the order 10 4 - 10 5 , which is quite high 
and might create numerical difficulties resulting in inaccurate solutions, unless the geometric quantities 
are approximated consistently. 

The vortex shedding is triggered by an asymmetric perturbation which consists of rotating the 
cylinder a short period of time in the clockwise and then in the counter-clockwise directions (ref. 5). 
The computation is performed within a limited non-dimensional time ( t < 40) in order to minimize the 
wall effects. It was found that a small amount of fourth-order artificial diffusion should be added for 
Re > 200 in order to smooth out the solution in the far wake region, where the grid is coarse. 
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The time evolution of the lift and the drag coefficients as a function of the Reynolds number is 
given in figure 19. In each case, following a relatively short transient flow, the onset of periodic flow 
dominated by a single frequency can be observed. The amplitude of the oscillations and the Strouhal 
number of the vortex shedding increase with the Reynolds number. Figure 20 compares the Strouhal 
number obtained in the present study with other numerical and experimental results (refs. 10, 19, and 30). 
Considering the complexity of the problem and the wide spread in the results, reasonable agreement is 
obtained, although the present results predict a slightly high Strouhal number for Re > 500. 

Figure 21 shows a partial view of the grid in the vicinity of the circular cylinder for the case 
of Re = 200. Figures 22(a) and 22(b) compare the instantaneous streamlines for five equally spaced 
instances along one vortex-shedding cycle with the results obtained by Rogers (ref. 31) for an identical 
case. The time is normalized by the period and t = 0 corresponds to the instant when the lift coefficient 
is maximal. The lift coefficient is minimal and the drag coafficient is maximal at t = j, while at 
t = ^ and t = | the lift coefficient vanishes and the drag-coefficient is minimal. The corresponding 
pressure fields for both solutions are shown in figure 23 for the same instances. The centers of the 
concentric pressure contours in the wake region correspond to a local minimum in the pressure. 

The periodic characteristics of the flow can be clearly seen. The flow fields at the beginning of the 
cycle ( t = 0) and at the end of the cycle (f = 1 ) are essentially identical. Moreover, the flow field at 
t = j is a mirror image of the flow field at t = 0. A similar relationship is found between the flow 
fields at t = ^ and t = when the instantaneous lift coefficient vanishes. A large separation bubble is 
created at the upper side of the cylinder at t = 0. The bubble induces a low pressure field, which causes 
a large positive lift coefficient, as well as maximal drag. The bubble ultimately separates from the body 
and is washed downstream ( t — ^). The lift coefficient is decreased to zero and the drag coefficient 
is minimal at that time, (fig. 19). At t = j a large separation bubble is created at the lower side, 
reversing the direction of the lift coefficient. This bubble is washed downstream as well and creates the 
lower vortex of the von Karman vortex street. The center of the shed vortices is characterized by a low 
pressure region. 

Rogers (ref. 3 1 ) solved an identical case using the same grid and time step but a different numerical 
method. His method is based on the artificial compressibility procedure with a fifth-order upwind 
scheme, second-order temporal accuracy and characteristic boundary conditions (refs. 6 and 7). Despite 
the differences in the numerical algorithms, good agreement is obtained, as shown in figures 22 and 23. 
The results compare particularly well in the regions near the circular cylinder, where the resolution of 
the grid is adequate in both solutions. The apparent discrepancy in the stream function contours at the 
center of the vortices is due to a small difference in the time of the corresponding plots (although the 
time increment is equal, the starting times of the shedding cycles are not the same since the triggering 
mechanism of the vortex shedding is different). The agreement deteriorates somewhat at the far wake 
region, where the mesh size is of the order of the cylinder’s diameter and truncation errors affect the 
solution differently for each case. 

Some dissipation of the downstream vortices can be observed in both solutions. The amplitude of 
the wake decreases and the steep pressure gradient in the vortices flattens out with increasing distance 
from the cylinder. This dissipation is mainly a result of the coarse grid at the far wake region (fig. 21). 
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Table 1 compares the Strouhal number and the drag and lift coefficients for the Re = 200 case 
with several other experimental and computational results (refs. 5, 10, 19, and 31-33). Gerrard’s results 
(ref. 32) summarize experimental results obtained from several sources. The scatter in both the numerical 
and the experimental results is quite large, demonstrating the complexity of the problem. The present 
computations are within the range of other available data. Note in particular the good agreement with 
the recent solution of Rogers (ref. 31). 

Figure 24 shows the effect of the Reynolds number on the instantaneous streamline contours for a 
time which corresponds to the maximal lift-coefficient at each Reynolds number. The laminar solution 
is computed in each case, although in the real flow the far wake does not remain laminar for Reynolds 
number higher than a few hundreds. The extent of the diffusion effects decrease with the increase of 
the Reynolds number. The attached separation bubble, where the diffusion effects are prominent, is 
the largest for Re = 100 and gradually decreases as the Reynolds number increases. The center of 
the attached separation bubble moves toward the lower side of the cylinder as the Reynolds number 
increases along with a downward bending of the streamlines. The wavelength of the wake decreases 
with increasing Reynolds number, in accordance with the increase in the Strouhal number (see fig. 20). 


Flow in a Square Duct with 90° Bend 

A square cross section duct with a 90° bend is solved as a test case for a three-dimensional internal 
flow. Flows in curved ducts occur in a wide range of practical applications such as aircraft intakes, 
turbomachinery blade passages, diffusers, and heat exchangers. A distinguishing characteristic of the 
flow in ducts with strong curvature is the generation of streamwise vorticity caused by the centrifugal 
forces which generate substantial secondary flow and redistribution of the longitudinal velocity. 

The geometry and the grid used in the present study are shown in figure 25. The symmetry of the 
problem across the plane x = 0.5 is utilized to save mesh points. The length of the side of the square 
cross section is set to one unit. The straight inflow and outflow sections before and after the bend are 
five units long, and the radius of curvature of the inner wall in the 90° bend is 1.8 units. Three different 
grids with mesh sizes of 21 x 11 x41,29x 15x51, and 41 x 21 x 61 were used for solving the problem. 
The finest mesh is shown in figure 25 along with an enlarged view of the grid at the crossflow plane. 

The Reynolds number is 790, based on the average velocity at the inlet. The fully developed 
laminar flow in a straight square duct is specified at the upstream boundary, while zero gradient of the 
velocity is given at the downstream boundary. The problem is solved time accurately, but only the 
steady state solution is presented here, for comparison with the experimental data of Humphrey et al. 
(ref. 34). The experimental data are limited to the distribution of the streamwise velocity over several 
sections of the duct. 

The computed streamwise velocity profiles over three different meshes are shown in figure 26. The 
streamwise velocity profile is presented in this figure for six cross sections along the duct. The location 
of these cross flow planes is shown in figure 25. At each plane, the streamwise velocity is given along 
two lines, x — 0.25 (fig. 26(a)) and x = 0.50 (fig. 26(b)). In most of the flow regions shown, the results 
of the two finer meshes are in very good agreement with the experimental results of Humphrey et al. 
(ref. 34), indicating that the solution is essentially grid-independent. 
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Figure 27 compares the streamwise velocity obtained from the finest grid computation with the 
numerical results of Rogers (ref. 3 1 ) (who used the same grid) and the experimental results of Humphrey 
et al. (ref. 34) at the same six crossflow planes along the duct. The agreement between the two 
numerical results is very good in most of the regions. The agreement with the experimental results is 
good, especially at the first four cross sections. Some discrepancy is found between the numerical and 
the experimental results at the two downstream planes. Nevertheless, the large peak of the streamwise 
velocity near the outer boundary is well predicted in the numerical computations. 

Figure 28 compares the present streamwise velocity solution with the numerical computation of 
Rogers (ref. 31) and the experimental results of Humphrey et al. (ref. 34) for the same cross sections. 
For each cross section, the distribution of the velocity along the x axis is given for five different radial 
locations. The agreement in the numerical results is remarkable, bearing in mind that totally different 
solution methods have been used. The agreement with the experimental results is also good in most of 
the flow regions. 

Figure 29 plots the crossflow velocity at the plane 9 — 90°. At that plane, a very complicated 
pattern of the secondary flow is already established. Three different vortices can be seen. A large and 
relatively strong vortex is found near the inner wall. A second vortex is found near the side wall, while 
a relatively weak vortex exits near the symmetry plane, at about the center of the duct. Figure 30 gives 
the pressure contours at several crossflow planes in the bend region and in the region downstream of the 
bend. The main vortex is generated in the bend and subsequently diffused in the downstream section. 

6. SUMMARY 

A method for solving the three-dimensional, unsteady, incompressible Navier-Stokes equations 
is presented for generalized curvilinear coordinate systems. Accuracy is achieved by finite-volume 
discretization on a staggered mesh along with consistent approximation of the geometric quantities. The 
formulation with the volume fluxes as dependent variables results in a simple extension of the staggered 
mesh approach to generalized coordinate systems and facilitates in satisfying the mass conservation 
equation to machine accuracy for all the computational cells. The fractional step method is combined 
with an approximate factorization of the momentum equations. Each step of the method is consistent 
and therefore the physical boundary conditions can be used. An efficient Poisson solver has been 
developed for generalized nonorthogonal coordinate systems using a consistent discrete approximation. 
The discrete approximation of the governing equations is second-order accurate in both time and space. 
However, in the present study the pressure computation is replaced by an approximate equation, reducing 
the temporal accuracy of the pressure to first order. 

Several two- and three-dimensional laminar validation cases have been presented for both steady and 
unsteady flows. Good agreement with other numerical and experimental results is obtained over a wide 
range of Reynolds numbers. Future work will include solution of time -dependent, three-dimensional 
flow fields for more complicated geometries. 
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APPENDIX A 


This appendix elaborates the calculation of S l that appear in the [/^-momentum equation at points 
different than where they are stored. All the formulas are based on linear averaging to satisfy the 
geometric constraint ( 7 ). Note that the same equations are used in the computation of S / and U as 
well. The terms that appear in the other momentum equations can be derived by cyclic permutation. 
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APPENDIX B 


This appendix elaborates the calculation of the pressure P at points other than the center of the 
computational cells. All the formulas are based on linear averaging and therefore the results are second 
order accurate. Only the terms that appear in the (/^-momentum equation are given. The terms 
that appear in the other momentum equations can be derived by cyclic permutation. Note that these 
equations can also be used for the computation of the volume (the inverse of the Jacobian) and the 
effective kinematic viscosity. The indices (z, j, k ) are omitted whenever possible. 

p i±r l ^p+Pi±>) 

p M = \(p+p^) 

p * ± { = 5< p + p *±i) < BI > 

P i±\,j±\ = 4^ + P i± 1 + P j± 1 + P i± U±l) 

P i±{,k±\ = ^( P + P i± 1 +P k± 1 +-^±1,^1) 
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APPENDIX C 


This appendix gives the operator R ^ for the momentum equation. 
by cyclic permutation. 
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APPENDIX D 


This appendix gives the implicit part of the diffusion terms of the momentum equation, D^(U^). 
The terms that appear in the other momentum equations can be derived by cyclic permutation. The 
operator is split into 

D t = D U + D Z* + D U < D1 ) 

Each sub-operator of j corresponds to a step of the approximate factorization method. 
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APPENDIX E 


This appendix gives the implicit part of the convection terms of the momentum equation, 
The terms that appear in the other momentum equations can be derived by cyclic permutation. 
The operator is split into 


c t - C U + c t,v + C U 

Each sub-operator of C^ i corresponds to a step of the approximate factorization method 
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APPENDIX F 


In the present implementation of the fractional step solution method with the explicit approximation 
of the convection terms, the computed velocity field is the exact solution of the discretized governing 
equations, provided the pressure field at the two previous time levels is the correct one and the system 
(30a) is solved exactly (without splitting errors). Yet, the computed pressure field is not the exact 
solution of the original equations, but only a first-order-accurate approximation in time. The following 
analysis does not apply to the case with 9 C > 0 . 

For ease of presentation, the equations are discretized only in time and written in a vector form 
equivalent to equations (18a) and ( 11 ) 

„n+l 

— — = G + 0 d V 2 u n+1 + 0 r V(AP) (FI) 


V • u n+ = 0 (F2) 

where G holds all the terms from the previous time levels n,n — 1 and A P = P n+1 — P n . Without 
losing generality, it is also assumed that a u = a p = 0. The velocity vector can always be split into 

u n+1 = u u} +u (j) (F3) 

such that U 0 is irrotational (see (F5)). Equations (FI) and (F2) can be split and solved exactly (since 
(FI) is linear) by the following fractional steps 

^=G + «jVV (F4> 

^ = e r V4> (F5) 

V • U0 = -V • Iiw (F 6 ) 

0 = 0 d V 2 u^ + 0 r V(AP - 0) (F7) 

with the proper boundary conditions. This splitting of the equations introduces a new unknown <p. Note 
that (F4) is an uncoupled equation, and (F5) and (F 6 ) should be solved as a coupled system for the 
unknowns and 0. Equation (F7) is a stand-alone equation for the unknown AP. 

In the present application of the fractional step method, only the steps (F4)-(F6) are solved. Instead 
of solving (F7) for the pressure, the first-order time-accurate substitution AP = 4> has been used. 
However, it is easy to verify that steps (F4)-(F6) yield an exact solution of the velocity, provided 
P n and P n_ * are known. 


In the practical application of the method, there is an error in the velocity field as well due to the 
fact that P n is not the correct solution for n > 1. However, the error in the velocity is still second-order 
in time (see (F5) and (F7)). 


PRECSDilYS r,*:.' 


r ruAiP 


i'AUA . ■ 




APPENDIX G 


The Laplacian-like operator L is assembled from the contributions: 


L — L ^ + Lrj + L ^ 
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The present appendix gives the L ^ term. The other terms may be obtained by cyclic permutation. 
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Table 1. Comparison of the lift ( C i) and drag (Cq) coefficients and the Strouhal number (St) for the 
flow over a circular cylinder at Re = 200 



Cd 

C L 

St 

Gerrard (ref. 32) (exp.) 

232 


0.18-0.20 

Wille (ref. 33) (exp.) 

1.3 



Lecointe and Piquet (ref. 5) 




2nd order 

1.46 ± 0.04 

±0.70 

0.23 

4th order 

1.58 ± 0.0035 

±0.50 

0.19 

Gresho et al. (ref. 10) 

1.76 ± 0.09 

±1.05 

0.21 

Braza et al. (ref. 19) 


±0.77 

0.20 

Rogers (ref. 31) 

1.33 ± 0.05 

±0.68 

0.19 

Present 

1.31 ± 0.04 

±0.65 

0.20 
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Figure 1. Definition of the primary cell. 



Figure 2. The computational cell (shaded) of the ^-momentum equation. 
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Figure 3 . The labeling of the points in the (77, £) plane for the Poisson equation solver. 






















Figure 5. Comparison of the streamlines, square Figure 6. Comparison of the vorticity contours, 
cavity problem. Re = 1000. square cavity problem. Re = 1000. 
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(b) 

Figure 7. Comparison of the velocity distribution (a) horizontal center line, (b) vertical center line for 
the square cavity case, Re = 1000. 
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Figure 11. Velocity direction for the polar cavity Figure 12. Comparison of the velocity distribution 
flow at Re = 350. along the radial lines 9 — — 20°, 0° and 20°. 
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(a) Re = 100 



(b) Re = 1000 


Figure 13. Comparison of the horizontal velocity component along the vertical centerline of the cubic 
cavity. 
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Figure 17. Time evolution of the separation length behind the circular cylinder at Re = 40. 
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Figure 18. Time evolution of the separation length behind the circular cylinder at Re = 550. 
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Figure 19. Time evolution of the force-coefficients on the circular cylinder at various Reynolds numbers. 



Figure 20. The dependence of the vortex-shedding Strouhal number on the Reynolds number for a 
circular cylinder. 
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